Astronomy & Astrophysics manuscript no. 10982man 


© ESO 2009 


July 9, 2009 





Time-dependent radiative transfer with PHOENIX 



D. Jack 1 , P. H. Hauschildt 1 , and E. Baron 



1,2 



On 
O 
O 
(N 



43 

6 

C3 



1 Hamburger Sternwarte, Gojenbergsweg 112, 21029 Hamburg, Germany 
e-mail: djack@hs.uni-hamburg.de; yeti@hs.uni-hamburg.de 

2 Homer L. Dodge Department of Physics and Astronomy, University of Oklahoma, 440 W Brooks, Rm 100, Norman, OK 73019- 
2061 USA 

e-mail: baron@ou.edu 



Received 8 September 2008 / Accepted 3 April 2009 



ABSTRACT 



^1- 

o 

OS 

o 



X 



Aims. We present first results and tests of a time-dependent extension to the general purpose model atmosphere code PHOENIX. We 
aim to produce light curves and spectra of hydro models for all types of supernovae. 

Methods. We extend our model atmosphere code PHOENIX to solve time-dependent non-grey, NLTE, radiative transfer in a special 
relativistic framework. A simple hydrodynamics solver was implemented to keep track of the energy conservation of the atmosphere 
during free expansion. 

Results. The correct operation of the new additions to PHOENIX were verified in test calculations. 

Conclusions. We have shown the correct operation of our extension to time-dependent radiative transfer and will be able to calculate 
supernova light curves and spectra in future work. 

Key words, radiative transfer - supernovae: general 



i— i 1. Introduction 



All types of supernovae are important for the role that they play 
in understanding stellar evolution and galactic nucleosynthesis 
and as cosmological probes. Type la supernovae are of particular 
cosmological interest, e.g., be cause the dark energy was discov- 
ered w ith Type la supernovae (Rie ss et al.ll9 98; Perlmutter et al.l 
Il999h . 

In dark energy studies, the goal now is to characterize the 
nature of the dark energy as a function of redshift. While there 
are other probes that will be used (gravitational lensing, baryon 
acoustic oscillations), a JDEM or Euclid mission will likely con- 
sider supernovae in some form. In planning for future dark en- 
ergy studies both from space and from the ground, it is im- 
portant to know whether the mission will require spectroscopy 
of modest resolution, or whether pure imaging or grism spec- 
troscopy will be adequate. Several purel y spectral ind i cators 



of peak luminosity have been proposed (iNugent et al 
Hachinger et al. 2006; Bongard et al. 2006; Bronder et 



tl. 
al. 



1995 



2008 



Foley et al. 20081; iLe Du et al.ll2QQ9l) 7 What is required is an em 



pirical and theoretic al comparison of both light curve shape lu- 
minosity indicators dPskovskiill 19771 iPhillipsI Tl 993b iRiess et all 
1996; Goldha beretal.1120011) and spectral indicators. 

To make this comparison one needs to know more about the 
physics going on in a supernova explosion and to be able to cal- 
culate light curves and spectra self-consistently. Thus, we need 
to extend our code to time-dependent problems. While our pri- 
mary focus is on Type la supernovae, the time-dependent radia- 
tive transfer code is applicable to all types of supernovae, as well 
as to other objects, e.g., stellar pulsations. 

In the following we present the methods we used to imple- 
ment time dependence. First we focus on solving the energy 
equation (first law) and then turn our attention to time-dependent 
radiative transfer. 



2. Equation of energy conservation 

To compute light curves we extended PHOENIX with a time- 
dependent solution of the non-grey, NLTE radiative transfer 
problem in special relativistic environments and a simple hy- 
drodynamic code to solve the free expanding case relevant to 
supernova light curves. The idea is to keep track of the conser- 
vation of energy for the gas and radiation together and to allow 
for different time scales for the gas and the radiation. 

For the time-dependent radiative transfer problem we 
extended our exist i ng r a diative transfer code PHOENIX 
(lHauschildt & Baronl 1 19991 l2004bl) . The change in the en- 
er gy density of a radiating m aterial is given by equation (96.15) 
in lMihalas & Mihalasl(ll984 

D d 
Dt dM r 

where E is the total energy density. All quantities are considered 
in the comoving frame. P r is not the pressure, but rather mechan- 
ical power on the sphere of a radius r. Equation Q] is only valid 
to first order in v/c, and thus lacks the full special relativistic 
accuracy of PHOENIX. This is adequate for the velocities found 
in supernovae. The total energy density of a radiating fluid con- 
sists of the sum of the energy density of the material, the energy 
density of the radiation field, the kinetic energy density of the 
material, and the gravitational energy density: 



E — En as + + Efcfa + Enyay. 

P 



(2) 



For supernovae in the free expansion phase, the gravitational 
energy density E grav is negligible since the potential is small 
in absolute value with the standard choice of zero at infinity. 
Homologous expansion is a reasonably good assumption for su- 
pernovae. The energy release by the decay of 56 Ni can influ - 
ence the dynamics of the expansion (Pinto & Eastman 2000). 
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IWooslev et al.l (120071) compared a study following this energy 
release to the results from assuming homologous expansion. 
Figure 2 in their paper shows the deviation and density variations 
can be as large as 10%. However, this is probably an upper limit 
due to the simple burning parameterization used in that study. 
Ultimately, when the deflagration to detonation transition is un- 
derstood, it will be important to revisit this issue and replace 1-D 
calculations with full 3-D calculations, which include the effects 
of clumps, as well as nickel bubble expansion. For now the ac- 
curacy of homologous expansion should be adequate, given the 
other uncertainties in the problem. 

With the assumption of homology, the velocity of a given 
matter element is then constant as is the kinetic energy density. 



Thus, we can neglect the kinetic energy term 



DE k , 
Dl 



So for our 



approach, we only have to consider the energy densities of the 
radiation field and the material. For the material, this includes 
effects such as an energy deposition due to radioactive decay of 
nickel and cobalt in an SN la envelope. 

The other possible cause of a change in t he energy density 
is the structure term. This term is given by dCooperstein et alj 
1986) 



^ (Pr + Lr) = ^ 4nr2[uiP + P ° )+F0] } 



8M r 



(3) 



where p is the pressure of the material and Pq the radiation pres- 
sure, u the velocity of the expanding gas, the radiative flux is 
represented by Fq, and the mass inside of the radius r of a layer 
is given by M r . The radiation pressure is a result of the solution 
of the detailed radiative transfer equation and given by 



An 

Po = —K, 

c 



(4) 



with K the second moment of the radiation field. The change of 
the energy density is given by the two quantities 



L r = 4nr F u 



and 



P r = Anr 2 u (p + P ) . 



(5) 



(6) 



If the atmosphere is in radiative equilibrium, the structure term 
is zero and the energy density stays constant if there is no addi- 
tional energy source and the atmosphere is not expanding. 

All the quantities required for the structure term can be de- 
rived from thermodynamics or the solution of the radiative trans- 
fer problem. We need the energy density of the material E gas 
and the energy density of the radiation field For the latter, 
we have to solve the radiative transfer equation for the radiation 
field and the radiative energy. We use our radiative transfer code 
PHOENIX to solve the time-independent radiative transfer equa- 
tion. The energy of the radiation field is given by 



E = 

c 

where J is the mean intensity and c the speed of light. 
The energy density of the material is given by 

3 p 3 R 

^gas ~ - 1 i 

2 p 2 m u 



(7) 



(8) 



with the mean molecular weigh m u and the universal gas con- 
stant R. The gas pressure is represented by p and the density by 



p. T stands for the temperature of the gas. The sum of the radia- 
tion and material energy density is then the total energy density 



Etotal — Eg as + 



Eq 
P ' 



(9) 



The change in this total energy density is given by Eq. 1 . So 
the equation to calculate the new energy density E new is given by 



d 

E„ ew = EoM - — — (L r + P r ) At + sAt. 
dM r 



(10) 



We now have all the needed equations to calculate a simple 
light curve. One problem for the calculation is that we can only 
determine the change in the total energy density for the next time 
step. However, the total energy change is divided into a change 
in the gas energy density and the energy density of the radiation 
field. To obtain the correct distribution of the gas and radiative 
energy, one has to iterate for each time step by solving the radia- 
tive transfer equation to compute the correct new temperature at 
the next time step. 

To get the correct new temperature we apply the following 
iteration scheme. The error in the energy density, E err is given 
by 



E current E target 



(ID 



^target 



Here, E target is the desired new total energy density, which is 
known from Eq. [TUl and E current is the total energy density ob- 
tained by Eq.|9]with the current temperature guess and the radia- 
tive transfer solution. Tests have shown that the error is almost 
linearly proportional to the temperature T. Therefore, a new tem- 
perature guess can be calculated for the next iteration step. The 
new temperature guess T„ eH , is obtained by 



Tnew — 



EerrPold E errojA T a 
Eerr Eerr a i.i 



(12) 



where T cur and E err are the current temperature guess and energy 
error. The variables T„id and E errM are the temperature and en- 
ergy error of the temperature iteration step before. With the new 
temperature guess we solve the radiative transfer equation again 
and check whether the total energy density is the desired one. It 
takes approximately four or five iteration steps to get the correct 
new temperature for a typical time step. The energy density is 
correct within an accuracy of 10~ 5 . 

3. Test light curves 

For the first test calculations we have solved the radiative trans- 
fer equation for a grey test atmosphere. Test model atmospheres 
are divided in 100 layers and we consider here time indepen- 
dent radiative transfer (in §[4] below we discuss time dependent 
radiative transfer). 

As a first test we applied our time evolution code to a static 
atmosphere. The test atmosphere is not expanding and no en- 
ergy sources are present. Inside the test atmosphere we used 
a "lightbulb" radiating with a constant luminosity to simulate 
the internal energy flow from a star. We assumed an approxi- 
mate temperature structure for t — and then let the atmosphere 
evolve towards radiative equilibrium. The atmosphere changed 
until it reaches steady state and the luminosity is constant in both 
space and time. The resulting final temperature structure should 
be identical to the structure for radiative equilibrium computed 
directly. 
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Fig. 1. Three different light curves for the evolution to a station- 
ary state. The three models have different luminosities produced 
by different inner "lightbulbs". 
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Fig. 2. The final temperature structures of the time evolution cal- 
culations. 



Fig. 3. The flux in every layer at each time step. 
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Fig. 4. The result for a sinusoidal varying lightbulb. Shown here 
is the luminosity in different layers. The phase shift between the 
lightbulb and the emergent flux is roughly n. 



In Fig. Q] we show the light curves for three different static 
models. We can observe the model on its way towards radiative 
equilibrium. All calculations were started from the same tem- 
perature structure. After a certain time, the radiative relaxation 
time scale, each atmosphere has the (constant) luminosity of the 
lightbulb throughout the configuration. 

The final temperature structure of the model atmosphere with 
the lightbulb with a luminosity of 10 31 erg/s is displayed in Fig. 
[2] Also shown is the result of a cal culation using the tempe rature 
correction procedure of PHOENIX dHauschildt et al.ll2.Q03l) to di- 
rectly compute the radiative equilibrium structure of the configu- 
ration. As one can clearly see, the resulting temperature structure 
of both methods agrees very well, and the maximum deviation 
of the temperature structure is less than one percent. Only the 
inner three layers have a deviation of up to 10 percent due to 
the implicit assumption of a diffusion approximation (and not a 
lightbulb) in the static calculation. 

The next test is to look at time varying atmospheres. As an 
example we considered an atmosphere with a sinusoidally vary- 
ing lightbulb inside. The resulting luminosity in each layer at 
each time point is plotted in Fig. [3] The luminosity in each layer 
is sinusoidal. It takes some time for the radiation to reach the 



outside boundary of the atmosphere and this results in a phase 
shift compared to the lightbulb. 

We now take a look at the luminosity in different layers. 
These are plotted in Fig. [4] One can see the phase shift of the 
sine in each layer. As expected, the luminosity from the central 
source needs time to move through the atmosphere. 

One can also see that the amplitude of the sine decreases with 
increasing radius. What one would expect is that the amplitude 
is the same in every layer. The radiation is moving outwards and 
the incoming varying luminosity from the lightbulb should move 
through every layer and therefore the amplitude is supposed to 
be the same in each layer. If we integrate the luminosity over a 
whole sine, it stays constant in every layer. In the plot one can 
see that the mean luminosity is at the same level, so the lumi- 
nosity is preserved. But why does the amplitude decrease? A 
possible explanation is that the radiation moving through the at- 
mosphere is going backwards when the sine is declining, so this 
smears out the amplitude. If the sine varying radiation is moving 
through a very thick atmosphere, the amplitude at the top of the 
atmosphere is finally flat. An observer sees the mean luminosity. 

For the next test we consider an atmosphere with an internal 
energy source. The initial structure is the radiative equilibrium 
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Fig. 5. The light curve of an lightbulb with an additional energy 
source. An energy source in each layer causes an increasing lu- 
minosity of the model atmosphere. 



structure of the static model with the lightbulb with a luminosity 
of 10 31 erg/s. We assume a constant energy deposition rate in 
each layer of the model atmosphere. The luminosity is expected 
to increase over time and towards the outside. Figure [5] shows 
a plot of the light curve of this test atmosphere. The luminosity 
increases in time because of the energy deposition 



3.1. Expansion 

In the case of SN la, homologous expansion is a good assump- 
tion after the initial breakout. In our test model each layer has a 
constant velocity, u oc r, to simulate a freely expanding envelope. 
The new radius r new of a layer is, therefore, determined by 



u ■ At + r \d 



(13) 



for a time step Af, while the layer is expanding with a velocity 
u. The radius before the new time step is r id- With the same 
assumption of homologous expansion the new density p new of a 
layer after at the new time step is determined by 



: Potd ■ 



r D id 



(14) 



where p \d is the old density. With the new radius and density, 
we are now able to calculate the new thermodynamics of the 
atmosphere. We assume the expansion of the supernova is an 
adiabatic process. The internal energy of the material changes 
due to work dW done during the expansion 



dE = dW = -pdV, 



(15) 



where p is the pressure and dV a volume change. The energy 
conservation equation considers the energy density. The change 
in the energy density for the adiabatic expansion is given by 

dW p 

dE adia = = -i-dV. (16) 

m m 

The mass m of a layer is given by its volume and density 

m = V-p. (17) 



Fig. 6. Light curve of an atmosphere that is just expanding. 



In our approach to a homologous expanding supernova of the 
type la, we consider each layer has a constant mass, so the 
derivative dV is given by 



m 

dV = -- T dp. 

P 2 

Together with the equation of state 

P = ~Tp, 
P 

we obtain the work for the adiabatic expansion as 

dW _ R T dp 
m P P 

For differences in a time step Af, the work is 



W, 



ad in 



w r , p 2 
A— = -Tin — 
m p \pi 



(18) 



(19) 



(20) 



(21) 



For the first test of a light curve for a supernova with an expand- 
ing atmosphere, we neglect the interaction between the layers, 
meaning the structure term is equal to the work of the adiabatic 
expansion so that 



D D (E 



(22) 



As we now consider expanding atmospheres, we use more 
supernova-like structures for the tests. We set the maximum ve- 
locity to 30000 km/s and use a radius of 5 x 10 15 cm. With this 
setup, we calculated a light curve for the expanding atmosphere. 

Figure [6] shows a plot of the light curve of supernova 
test atmosphere that is simply expanding. We considered time- 
independent radiative transfer and a grey atmosphere. As one 
can see, the luminosity is decreasing, because the atmosphere is 
cooling down adiabatically. 

Now we test a setup more closely resembling a real super- 
nova light curve. Therefore, we take an initial atmosphere struc- 
ture and add an energy source (radioactive decay) in each layer. 
The energy source exponentially decreases to simulate declin- 
ing activity of the radioactive species. Figure Q shows the plot 
of the light curve of this test, resulting in a light curve with a 
supernova-like shape. We have a rising part of the light curve at 
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Fig. 7. Light curve of an atmosphere that is expanding and has 
an energy source. It has the typical shape of a light curve of a 
SN la. The luminosity rises due to the energy deposition. After 
the maximum we see the decline resulting from the expansion. 

the beginning because of the energy deposition. After the max- 
imum, the luminosity decreases due to the ongoing expansion 
and decreasing energy deposition. Of course this light curve is 
far from correct because the assumption of a grey atmosphere is 
a bad assumption for an SN la. But the tests show that the code 
behaves as expected. 

3.2. Entropy 

We calculated the entropy to test the code. Because energy 
is moving through the atmosphere, entropy is not conserved. 
Nevertheless, testing entropy conservation is an important test 
of the correctness of the code. Therefore, we calculated the en- 
tropy change in the case of pure adiabatic expansion. For testing, 
we neglected interaction between layers and gamma ray deposi- 
tion and calculated the entropy change for a time step for the just 
expanding case. 

To be consistent with our hydrodynamics equations, we de- 
duced the entropy from the first law of thermodynamics. A 
change in the entropy during a time step is therefore given by 



AS 



mR 2 n \T\j 2\fx 2 



1 



-Ilnfe 



(23) 



where 2 is the index of quantities at the new time and 1 the one 
of the old. For the integration of the temperature and the density 
term, we kept fi fixed. This is not correct, but it is simpler to 
solve and the resulting differences for the entropy are small. 

The entropy given here is only the entropy of the gas, so to 
check that the expansion worked right we had to neglect the ra- 
diation energy density in our energy density conservation equa- 
tion. Furthermore we ignored the interaction between the layers 
and neglected the structure term. We now have only a change in 
the energy due to the adiabatic expansion. 

With this setup we calculated a few time steps and checked 
the entropy. Even for a long time step of 1000s the entropy stays 
almost constant. The relative change of the entropy was » 10" 6 . 

4. Time-dependent radiative transfer 

So far, our radiative transfer code only solves the time- 
independent radiative transfer equation. For an implementa- 



tion of the time dependence in the radiative transfer itself, the 
spherical symmetric special relativistic radi ative transfer equa- 
tion ( SSRTE) for expanding atmospheres (Hau schildt & B aron 
1999) is extended so that the additional time dependent term is 
given by 



c Ot 



(24) 



where /3 — j- is the velocity in units of the speed of light c and 
y — (1 - f3 2 y l l 2 is the usual Lorentz factor. Here, I is the in- 
tensity, n the cosine of the angle between the radial direction 
and the propagation vector of the light. Using the notation of 
Hauschildt & Baron (2004a), the comoving frame SSRTE with 
the additional time dependent term is given by 



dl dl dl dAI 

a '"37 + °rT + a ^ + a -i"3T + 4fl <i / = *7 ~X l 
dt Or OjJ. OA 



(25) 



where 77 is the emissivity and \ the extinction coefficient. The 
wavelength is represented by A. The additional time dependent 
coefficient is given by 



Along the characteristics the equation has the form 
dl, dl dAI 

~r + ai T: + a '~^T = T ii-0a + 4a/)// 

ds ot OA 



(26) 



(27) 



where ds is a line element along a (curved) characteristic and //() 
is the specific intensity along the characteristic at point s > 0(s = 
denotes the beginning of the characteristic). For a de finition of 
the other coefficients see Hauschildt & Baron (2004a). 



4.1. Discretization of the time derivative 

We used the fi rst di scretization as described in 
lHauschildt & Baronl d2004a) and added the time dependent 
term. We discretized the time-dependent, as well as the wave- 
length, derivative in the SSRTE with an fully implicit method. 
The discretization of the time dependent term is given by 



ai_ 

at 



h -i t , 



(28) 



so the SSRTE including the time discretization can be written as 



dl Ail - Ai-ih,_ 

h a i 

ds A t - At-i 



■ a/—^ = m , - ( X a, + 4^)7, (29) 



where I is the intensity at wavelength point Ai and time point tj. 
We define the optical depth scale along the ray as 



dr = y + a A K + — — 1 + — = -Ms. 

X A \ A,-A^j At X 



A, - A h 

Introducing the source function S 
dr x\ X h ~ 



T]/x, we get 

= I-S. 



a, 1 



(30) 



(31) 



We now have a modification of the source function and a 
new definition of the optical depth scale along a ray. With 
this redefinition of the optical depth and the source function, 
one can now proceed with the formal solution as described in 
lHauschildt & Baronl d2004al) . 
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Fig. 8. Luminosity of each layer and time step of our test atmo- Fig. 10. The light curve of an atmosphere where a small pertur- 
sphere. bation moves outwards. The two plotted light curves of layer 110 

were calculated with different time steps. 
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Fig. 9. The light curves of a few layers. In this case the atmo- 
sphere has an additional light source switching on inside. The 
radiation needs time to get to the outer layers. 



4.2. Test calculations: time-dependent radiative transfer 

For the test of the time-dependent radiative transfer, we used a 
static atmosphere structure to see the direct effects of the time 
dependence of the radiative transfer. Therefore, the tempera- 
ture, radius, and density are all constant in time. For the test we 
then changed the inner boundary condition for the radiation (the 
"lightbulb") to initiate a perturbation of the radiation field, which 
then moves through the atmosphere via the time-dependent radi- 
ation transfer. 

For the first test we switched on an additional light source 
inside the atmosphere. This light source has a luminosity of 10 9 
times higher then the inner lightbulb. In Fig. [8] we show the re- 
sult of the time-dependent radiative transfer calculation. One can 
see that the light propagates outwards through the atmosphere. 
It takes a while before the additional light has propagated every- 
where throughout the atmosphere. In Fig.[9]the light curves of a 
few layers are shown. One can see that the radiation needs time 
to get to the outer layers. 

We compared our time scale to the radiation diffusion time 
scale. Assuming a random-walk process the mean free path for 



a photon is given by A p - I, where x 1S the mean opacity. For a 
travel distance I the time t„ a photon needs is given by 



3 C*' 



(32) 



where c is the speed of light. For the mean opacity x we use d 
the Rosseland mean. In our test case the atmosphere was divided 
into 128 layers. The mean opacity ;f ranges from 3 • 10~ 19 ctrT 1 
in the outer parts to 3 • 10~ 14 ctrT 1 in the inner parts of the at- 
mosphere. The distance / is the thickness of each layer, and it 
ranges between 6 • 10 12 cm and 3 • 10 13 cm. The result is that the 
diffusion time for a photon through the whole model atmosphere 
is about 3 • 10 3 s. As one can see in our plots our time scale is 
roughly 10 5 . The assumption of a diffusion through the atmo- 
sphere is only valid for optically thick regions. Another problem 
is the choice of the right mean opacity. Considering this our es- 
timate of the time scale is reasonable. 

An important test is also to check that the results of the time- 
dependent radiative transfer calculation depend on the size of 
the time step. We tested this with a model that has a small per- 
turbation inside that is moving outwards. This was calculated 
with two different time steps. In Fig. [I0]the results of calcula- 
tions with two different time steps are shown. As one can see, 
the result does not depend on the size of the time step. 

In the next test we put a sinusoidally varying light bulb inside 
of our test atmosphere. In Fig.Qj] we show a plot of the luminos- 
ity for each time step. After a while, the luminosity of the whole 
atmosphere varies sinusoidally and steady state is reached. The 
luminosities in different layers are shown in Fig. [12] One can see 
a phase shift of the sine curve, due to the time required for the 
radiation field to propagate through the atmosphere. 



5. Conclusion 

We have implemented a time evolution code into our general- 
purpose model atmosphere code PHOENIX, which keeps track of 
the energy conservation. Because a homologous expansion is a 
good assumption for supernovae in general and particularly for 
Type la supernovae, we considered adiabatic free expansion for 
our code. With first test calculations, we reproduced the expected 
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Fig. 11. Luminosity in each layer and time step. The light source 
inside the atmosphere varies sinusoidally. 
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dEastman & Pintoll993tlBaron et al.ll996l) . One challenge is the 
computation time needed for a whole light curve with more com- 
plex NLTE radiative transfer. 
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Fig. 12. Flux at different radii of a sinusoidally varying atmo- 
sphere. 



behavior of the test cases. Static atmospheres adopted the lumi- 
nosity of an internal lightbulb and perturbations of the lightbulb 
moved outwards in time. We also calculated a light curve that 
had the shape of a typical supernovae la light curve. 

To complete the physics of time-dependent processes in 
a supernova atmosphere, we extended our code even further. 
PHOENIX now solves the time dependent radiative transfer equa- 
tion. We presented the new discretization scheme of the time 
dependent term. We checked our code with a series of tests. A 
perturbation can be followed on its way through the atmosphere. 
We ran a test model with a sinusoidal source. In steady state the 
whole atmosphere was varying sinusoidally, responding to the 
source. A phase shift between the inner and outer layers could 
also be observed. All tests indicate that our extended code works 
fine. 

Future work is to calculate realistic light curves for Type la 
supernova hydro models. As a first step we can consider the 
atmosphere to be in LTE, but including detailed opacities and 
treating full line blanketing. 

Also with our time-dependent radiative transfer, we can 
address a longstanding debate about the importance of time 
dependence in calculating the spectra of Type la supernovae 



